library(maptools)
## Loading required package: sp
## Checking rgeos availability: TRUE
library(lattice)
library(spdep)
## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
## Loading required package: sf
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(sp)
library(rgdal)
## rgdal: version: 1.5-18, (SVN revision 1082)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: D:/Program Files/R-4.0.3/library/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: D:/Program Files/R-4.0.3/library/rgdal/proj
## Linking to sp version:1.4-4
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
library(tmap)
library(ggplot2)
library(gridExtra)
library(gstat)
library(OpenStreetMap)
library(spacetime)
library(knitr)
library(tmap)
library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --
## √ tibble  3.0.4     √ dplyr   1.0.2
## √ tidyr   1.1.2     √ stringr 1.4.0
## √ readr   1.4.0     √ forcats 0.5.0
## √ purrr   0.3.4
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::combine() masks gridExtra::combine()
## x dplyr::filter()  masks stats::filter()
## x dplyr::lag()     masks stats::lag()
library(spdep)
library(sf)
library(OpenStreetMap)
library(spatstat)
## Loading required package: spatstat.data
## Loading required package: nlme
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
## 
##     collapse
## Loading required package: rpart
## Registered S3 method overwritten by 'spatstat':
##   method     from
##   print.boxx cli
## 
## spatstat 1.64-1       (nickname: 'Help you I can, yes!') 
## For an introduction to spatstat, type 'beginner'
## 
## Note: spatstat version 1.64-1 is out of date by more than a year; we strongly recommend upgrading to the latest version.
## 
## Attaching package: 'spatstat'
## The following object is masked from 'package:gstat':
## 
##     idw
## The following object is masked from 'package:lattice':
## 
##     panel.histogram
library(here)
## here() starts at D:/CASA/STDM/practical 1
library(sp)
library(rgeos)
## rgeos version: 0.5-5, (SVN revision 640)
##  GEOS runtime version: 3.8.0-CAPI-1.13.1 
##  Linking to sp version: 1.4-4 
##  Polygon checking: TRUE
library(maptools)
library(GISTools)
## Loading required package: RColorBrewer
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:spatstat':
## 
##     area
## The following object is masked from 'package:dplyr':
## 
##     select
library(tmap)
library(sf)
library(geojson)
## 
## Attaching package: 'geojson'
## The following object is masked from 'package:graphics':
## 
##     polygon
library(geojsonio)
## 
## Attaching package: 'geojsonio'
## The following object is masked from 'package:base':
## 
##     pretty
library(tmaptools)

dt <- st_read("data/factorscore.geojson")
## Reading layer `factorscore' from data source `D:\CASA\STDM\practical 1\Data\factorscore.geojson' using driver `GeoJSON'
## Simple feature collection with 779 features and 5 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 523850.3 ymin: 177849.8 xmax: 531172.2 ymax: 183893.7
## geographic CRS: WGS 84
dt <- dt %>%
  st_set_crs(., 27700)  #adjust the proj
## Warning: st_crs<- : replacing crs does not reproject data; use st_transform for
## that
##dt <- st_transform(dt, 4326)


# the spatial matrix
W_cont_el <- poly2nb(dt, queen=T)

#Queen's Case

W_cont_el_mat <- nb2listw(W_cont_el, style="W", zero.policy=TRUE)




#calculate Getis-Ord Gi* scores

lg1 <- localG(dt$region_1, listw=W_cont_el_mat, zero.policy=T)





dt$lg1 <- lg1[]



dt <- as(dt, 'Spatial')
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Airy 1830 ellipsoid in CRS
## definition
## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum OSGB 1936 in CRS definition
tmaplondon <- dt %>%
  st_bbox(.) %>% 
  tmaptools::read_osm(., type = "osm", zoom = NULL)  ##london osm map
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
breaks1<-c(-1000,-2.58,-1.96,-1.65,1.65,1.96,2.58,1000)

GIColours<- rev(brewer.pal(8, "RdBu"))


tmap_mode("plot")
## tmap mode set to plotting
#now plot on an interactive map


tm <- qtm(tmaplondon) + tm_shape(dt) + 
  
  tm_polygons("lg1",
              style="fixed",
              breaks=breaks1,
              palette=GIColours,
              midpoint=NA,
              alpha = 0.7,
              title="Getis-Ord Gi* of factor0") + tm_layout(legend.position = c("right", "top"), legend.bg.color = "white", legend.bg.alpha = 0.7)+  
  tm_compass(position = c("left", "bottom"),type = "arrow") + 
  tm_scale_bar(position = c("left", "bottom"))

#tmap_save(tm, filename = "map1.png")
tm

###############################

lg2 <- localG(dt$region_2, listw=W_cont_el_mat, zero.policy=T)
dt$lg2 <- lg2[]
tm1 <- qtm(tmaplondon) + tm_shape(dt) + 
  
  tm_polygons("lg2",
              style="fixed",
              breaks=breaks1,
              palette=GIColours,
              midpoint=NA,
              alpha = 0.7,
              title="Getis-Ord Gi* of factor1") + tm_layout(legend.position = c("right", "top"), legend.bg.color = "white", legend.bg.alpha = 0.7)+  
  tm_compass(position = c("left", "bottom"),type = "arrow") + 
  tm_scale_bar(position = c("left", "bottom"))

tm1

#tmap_save(tm1, filename = "map2.png")

#########################

lg3 <- localG(dt$region_3, listw=W_cont_el_mat, zero.policy=T)
dt$lg3 <- lg3[]
tm2 <- qtm(tmaplondon) + tm_shape(dt) + 
  
  tm_polygons("lg3",
              style="fixed",
              breaks=breaks1,
              palette=GIColours,
              midpoint=NA,
              alpha = 0.7,
              title="Getis-Ord Gi* of factor2") + tm_layout(legend.position = c("right", "top"), legend.bg.color = "white", legend.bg.alpha = 0.7)+  
  tm_compass(position = c("left", "bottom"),type = "arrow") + 
  tm_scale_bar(position = c("left", "bottom"))

tm2

#tmap_save(tm2, filename = "map3.png")

############################



dt$Kmeans_Cluster = as.character(dt$Kmeans_Cluster)


tmc <- qtm(tmaplondon)+tm_shape(dt) + 
  tm_polygons("Kmeans_Cluster", alpha = 0.85, palette = 'Pastel1', title="Result of Cluster", border.alpha = 0.5
              )+ tm_layout(legend.position = c("right", "top"), legend.bg.color = "white", legend.bg.alpha = 0.7)+  
  tm_compass(position = c("left", "bottom"),type = "arrow") + 
  tm_scale_bar(position = c("left", "bottom"))


#tmap_save(tmc, filename = "mapc.png")

tmc